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Open systems dynamics: Simulating master equations in the computer 

Carlos Navarrete-Benllodlf] 

Abstract 

Master equations are probably the most fundamental equations for anyone working in quantum optics in the presence of 
dissipation. In this context it is then incredibly useful to have efficient ways of coding and simulating such equations in the 
computer, and in this notes I try to introduce in a comprehensive way how do I do so, focusing on Mat lab, but making it general 
enough so that it can be directly translated to any other language or software of choice. I inherited most of my methods from 
Juan Jose Garcia-Ripoll (whose numerical abilities I cannot praise enough), changing them here and there to accommodate 
them to the way my (fairly limited) numerical brain works, and to connect them as much as possible to how I understand 
the theory behind them. At present, the notes focus on how to code master equations and hnd their steady state, but I hope 
soon I will be able to update them with time evolution methods, including how to deal with time-dependent master equations. 
During the last 4 years I’ve tested these methods in various different contexts, including circuit quantum electrodynamics, the 
laser problem, optical parametric oscillators, and optomechanical systems. Comments and (constructive) criticism are greatly 
welcome, and will be properly credited and acknowledged. 
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I. ON THE STRUCTURE OF MASTER EQUATIONS AND STEADY STATES 


Let me start by briefly introducing in a greatly simplified manner the concepts of master equation and steady states. 
Consider a system endorsed with a Hilbert space H of dimension d (since we are interested in doing numerics, we 
will always assume that d is finite, what might need truncating the Hilbert space dimension when this is infinite in 
reality). We say that the system is open when it is part of a larger space with which it exchanges energy, information, 
etc..., and generically we call environment to the rest of this larger space. In many situations, most notably when 
the environment is much larger or evolves much faster than the system, it is possible to describe the dynamics of the 
latter via a linear differential equation for its individual state, which of course is generally mixed since eliminating 
the environment means loosing information, hence requiring a description in terms of a density operator p. We call 
master equation to the evolution equation for the system’s density operator, and in the following we will be guided 
by the generic form^ 


dp 

dt 



+ T{2Jpr -JUp- pjU) = C[p], 


( 1 ) 


where H is an Hermitian operator containing the system Hamiltonian and coherent or reversible exchange processes 
with the environment, while J is a so-called jump operator (with associated rate L > 0 ) which describes irreversible 
processes such as excitations which are lost to the large environment never to come back to the system. C is then an 
linear map usually called Liouvillian superoperator, whose name comes from the fact that it acts on operators to give 
operators. 

Consider a basis {\n)}n=i,2,...,d^ which allows us to represent the density operator as 


d 

P= Yi Pnm\n){m\, ( 2 ) 

nm=l 

with pnm = {n\p\m). The master equation, once projected into this basis, just provides a linear system of ordinary 
differential equations for the components of the density matrix, that is 


dp 


nm 


dt 


d d d 

^ {ipnkHkm - iHnkPkm) + Y (‘^^^nkPklJLl ” ^J^JuPl m ^ PnkJlkJlm) — nm^klPkl 5 

k=l kl=l kl=l 


( 3 ) 


with 


Lnm]kl — iHlm^kn “ 1 “ JnkJml J^nk^ml J^lm- 

In a more compact notation, it is customary to take the columns of the density matrix, and pile them one below the 
previous one, transforming the matrix into a vector 


p Col(pii, P21 5 • • • 5 Pdl 5 Pi2 5 P221 “‘1 Pd 2 5 • • • 5 Pld-) P 2 d') • • • 5 Pddj) 5 


( 4 ) 


and the multidimensional array {Lnm\ki}nmki=i,2,...,d iiifo a matrix L, so that the previous equation is turned into 
the linear system 


with solution 



( 5 ) 


p{t) = e’^*p(0). 


(6) 


It is interesting to note the correspondence between the elements of the density matrix, and the elements of its 
veetorized form: (p)n+(m-i)d = Pnm- In the next section we will see that these expressions are much more than just 
a convenient way of reordering things. 


^ We will consider a single jump operator for notational simplicity, but everything we’ll do is generalized straightforwardly to the general 
irresversible term Yj ^pJj ~ j]jjP ~ pj]jj)^ even to terms of different form. 
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From a practical point of view, if the Hilbert space dimension is not too large, then can be efficiently evaluated, 
and the main problem consists in how to write the matrix L in an easy way, starting from the expression of the 
Liouvillian in the master equation 0. This is where superspace enters into play, and we will explain in the next 
section how it allows for a simple way of coding L in the computer (Matlab in particular). 

Let us now pass to discuss the important concept of steady state. For problems without selective measurements 
involved in the system+environment, the master equation must map states into states; this means that it is a trace 
preserving differential map, so that the condition tr{p} = pn + p22 + ... + pdd = 1 is satisfied at all times, and 
hence the equations of <§ are not independent, but satisfy the constrain pn + P22 + ... + pdd = 0, the dot denoting 
time-derivative. This makes the rows (or columns) of the matrix L linearly dependent, what ensures det{L} = 0 and 
hence that it exists at least one eigenvector with zero eigenvalue, which we will denote by po, satisfying Lpo = 0, where 
0 is a vector of zeros. The corresponding operator po is called the steady state of the system, since in the absence of 
any other zero eigenvalue, this is the state towards which the system tends to as time evolves (the trace-preservation 
condition ensures also that all the other eigenvalues have negative real part, and hence for long times only po survives). 

One useful way of finding this steady state is as follows. First, given its defining equation Lpo = 0, where 0 is a 
vector of zeros, we replace one of the equations coming from the evolution equation of some some chosen diagonal 
element pu by the normalization condition 7(pii + P22 + ... + pdd) = 7 , which can be additionally multiplied by any 
number 7, what is sometimes useful for numerical purposes; this means replacing the row number / + (/ — l)d of L 
by a vector of 7’s in the elements multiplying the diagonal elements of po, obtaining a new matrix Lq, also replacing 
the 0 vector by a vector wq containing a single non-zero entry 7 at position I ^ {I — l)d. The steady state can then 
be found by solving the linear system Lopo = wq^ for example by inversion: po = Later we will learn how to 

do this explicitly in Matlab. 


II. THE MASTER EQUATION IN SUPERSPACE 

The space where the density matrix is turned into a vector and the Liouvillian into a matrix is usually called 
superspaee. Having operators as its elements, superspace can be defined formally as the tensor product of the Hilbert 
space and its dual, which indeed has vector space structure when endorsed with the trace product. As we will see in 
the next section, this gives us a simple way of representing superoperators by using simple tools of computer programs 
such as the Kronecker product, which is a built-in operation in both in Matlab and Mathematica. Instead of recalling 
the dual space, it is computationally more convenient to define superspace in a slightly simpler way: given an operator 
O with matrix elements {Onm}n,m=i,2,...,d5 we just associate to every index a fictitious d-dimensional Hilbert space 
with basis {\n)}n=i^2,...,d^ which we describe the operator as a vector 

d 

|0) = y; Onm\n) O \m). (7) 

nm=l 

We will say that \ 0 ) is the abstract superspace representation of the operator O. Note that its representation in the 
basis {|n) (8) |^)}n,m=i,2,...,d of superspace is the vector O obtained by piling up the columns of the matrix formed by 
the elements Onm one below the previous one, but only provided that the order of the superspace basis {|p)}p=i,2,...,d2 
is chosen as^ {|n + (m — l)d) = |n) (8) 

Consider now two operators A and H, and a superoperator S acting on an operator O as 5 [ 0 ] = AOB = 
^nmkl=l OnmA knBmi\k){l\j expression which reads in superspace as 

d 

I] OnmAknBml\k)<E>\l). (8) 

nmkl=l 


Taking into account that A and B are operators defined in the original d-dimensional Hilbert space H, so that they 
act on basis elements of the new fictitious spaces in the usual way (e.g., A\n) = ^kn\k))^ we can alternatively 


^ Convince yourself of this fact through some simple examples. For example, the first element of the second column, O12, should correspond 
to the d + 1 element in the superspace vector, ( 0 )(^_^i, and this is precisely what the map (n, m) n + {jn — l)d provides. 
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write 


d 

|5[0])= ^ O„™(i|n)0B^|m)) = (i0B^)|d), (9) 

nm=l 

showing that in superspace the action of operators on the left (right), corresponds to actions of the (transpose) 
operator on the first (second) fictitious Hilbert space. 

Hence, in superspace the master equation can be written as 

^ 1/5) = [-i{H ® /) + i(/ ® + 2r( j ® j*) - r( jt j ® /) - r(/ ® j*)] \p), (lo) 

where I is the identity operator. 


III. IMPLEMENTING THE SUPERSPACE IDEAS IN MATLAB 

Let’s pass now to discuss how to implement the previous ideas in one particular program, Matlab, although similar 
tricks can be used in Mathematica, for example. 

First, let us remind the notation, since I wouldn’t be surprised if everyone is lost on it at this point; to complicate 
things a bit, we will even need to introduce some more. The (column) vector representation of the basis elements 
{|n)}n=i,2,...,d of the Hilbert spaces H or T will be denoted by {vn}n=i, 2 ,...,d^ with self-representation elements 

(v^)^ = {m\n) = Smn- Similarly, the superspace basis {\p)}p=i^2,...,d‘^ will have a vector representation {vp}p=i^2,...,d‘^^ 
with self representation (vp)^ = Sqp. Given an operator O, its matrix elements are denoted by Onm = (n| 0 |m), and 
the d X d matrix that they form by O, which is nothing but the matrix representation of the operator in the chosen 
basis. The ‘vectorized’ form of this matrix, that is, the vector formed by piling up the columns of the matrix one 

below the next, is denoted by O, and, as explained above, it can be seen as the representation of the operator in 

superspace, which we will still denote as |0), provided that the basis elements of superspace are ordered in the proper 
way. Hence, summing up, O is the abstract notation for the operator acting on the original space H and \0) the 
one for the operator defined in superspace T ^ with corresponding matrix and vector representations O and O, 
respectively. As for superoperators, take the Liouvillian as an example, we will refer to them as C in calligraphic font^ 
when talking about them in an abstract way, and L in blackboard font when referring to their matrix representation 
in superspace. For example, the right hand side of master equation 0 reads C[p] in an abstract way, and as Lp once 
represented in superspace. 

Let’s start from the basics of coding things in Matlab. The matrix representation of an operator O can be written 
in Matlab (known their matrix elements Onm) as 

O = [Oil, Oi2, •••, Old] 021,022, •••, 02d] •••; 0dl,0d2, •••, Odd], (11) 

where the commas can be replaced by a space. If we want to work with sparse matrices to save memory (useful for 

Hilbert spaces with large dimension, e.g., d > 10 ), we can do so by replacing the matrix O by^ sparse(O); once in 
sparse form, we can always come back to the non-sparse one as full(O). 

Let’s talk about basic matrix operations. Element Onm is accessed as 0 (n,m), while the whole column m can be 
accessed as 0 (:,m), and similarly for row n, 0 (n,:). We can access its n-th diagonal as diag(0,n) which generates 
a column vector with the desired diagonal. Given another matrix Q, O zb Q is the matrix sum or difference, 0 *Q 
is the matrix multiplication, 0 .*Q is the element-by-element multiplication, 0 \Q is the matrix multiplication of 
and Q, and O/Q is the matrix multiplication of O and Q“^, where the inverse of O can also be obtained as 
inv(O) (but it is not recommended by Matlab, since it is slower than I/O or 0 \ 1 ). We can find the determinant 
and trace as det(O) and trace(O). The Hermitian conjugate of O is obtained in Matlab as O', while O.' generates 
the transpose of O. The exponential matrix is obtained as expm(O), while exp(O) just exponentiates the elements 
of the matrix individually. As for the eigensystem, eig(O) generates a vector with the eigenvalues of O, while if 
we write [V, D] =eig(0), the eigenvectors are codified as columns of V and D is a diagonal matrix containing the 


^ Not to confuse with the notation for Hilbert spaces, for which we also use calligraphic font, but it should be clear from the context when 
we mean one or the other. 

^ From now on, Matlab functions will be highlighted by using typewriter font. 
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corresponding eigenvalues. This operation cannot be used when O is sparse, in which case we need to use eigs(O), 
which by default gives the 6 eigenvalues with largest magnitude; if we want a different number of eigenvalues, say 
N, we can write eigs(0, where x = L (S) means that we want the eigenvalues with the largest (smallest) 

magnitude if ^ = M, or real part if ^ = R. It is very useful to type help F in Matlab’s command window to get more 
info about some function F (for example, try help eig to find what else can be done with eig). 

The vectorized form of the operator is obtained as O = 0(:). This is one of the reasons why we chose to pile columns 
instead of rows when vectorizing matrix representations of operators: in Matlab this operation is written with a single 
order, while in the case of rows, we would need to first transpose the matrix, and then give the order. On the other 
hand, given the matrix in vectorized form, we can always bring it back to matrix form as O = reshape(0, d, d). 

The matrix representation of the identity operator I can be written as I = eye(d), or I = speye(d) when working 
with sparse matrices. The basis vector can be defined as = I(^^)- Similarly, defining the identity matrix in 
dimensions I = eye(d^), the basis vector Vp in superspace is obtained as Vp = I(:,p). 

Let’s move on to the tensor product operation. It is customary in quantum mechanics to represent the tensor 
product of two operators or vectors as the Kronecker product of their representations. In particular, given the matrix 
representations A and B of two operators A and 5 , their Kronecker product is defined as 


/ AuB 

A12B • • 

• ^idB\ 

A21B 

A22B • • 

• A 2 dB 


Ad 2 B • • 

• ^dd'^ ) 


( 12 ) 


and this is usually the representation chosen for the tensor product operator A 0 B. However, it is important 
to understand that this is just one possible representation of the tensor product, corresponding to one particular 
ordering of the tensor product basis {|n) ( 8 ) \m)}n,m=i,2,...,d in the composite Hilbert space (superspace in our case). 
More concretly, note that according to the previous definition, given the vector representation of the basis element 
\n) G B (or J^), the Kronecker product representation of the basis element |n) (8) |m) G H (8) H (or B 0 B) generates a 
vector with a single nonzero entry at position m + (n — l)d, that is, the superspace basis vector^ However, 

as explained in the previous section, we would like to associate instead the superspace basis vector to 

the tensor product basis element |n) ( 8 ) |m), what means that we will not be using the usual Kronecker product 
representation of the tensor product, but one in reversed order: given the matrix representations A and B of two 
operators A and H, the representation of their tensor product A 0 B is taken as their Kronecker product in reversed 
order, that is. 


/BiiA 

B12A 

• 5 idA\ 

-B21A 

H22A • • 


V Ba-vA 


• Bdd-A j 


(13) 


With this choice, the representation of the basis element \n) ( 8 ) \m) corresponds to Vn-\-(m-i)d as we wanted to. 

The Kronecker product is already implemented in Matlab through the operation kron (and same in Mathemat- 
ica), which preserves the sparse character of the matrices. Hence, we can generate the vector representation of 

|n + (m — l)d) = \n) ( 8 ) |m), by applying the kron operation in the reversed order = kron(v^, v^). 

Consider now two operators A and H, and a superoperator S which acts on a third operator O as 5 [ 0 ] = AOB. 
As explained in the previous section, in superspace this is rewritten as | 5 [ 0 ]) = (A ( 8 ) B'^)\ 0 ) in an abstract way, 
expression which can be represented in the basis of superspace {\p)}p=i^2,...,d‘^ as SO = kron(B.', A)*0(:) in Matlab 
code. Hence, the matrix representation of S in superspace is written in Matlab as S = kron(B.', A). 

Let me remark that the choice of using the reversed kron order for the representation of the tensor product has 
been made for convenience in Matlab (to create the vectorized form of any operator with a single instruction, and 
for more things that will appear in the next section when dealing with composite Hilbert spaces). However, in other 
languages such as Mathematica, it can be better to stick to the traditional Kronecker product representation of the 


^ Convince yourself of this fact by working out some examples, e.g., in dimension 3 (d = 3), the Kronecker product of vi = col(l,0,0) 
and V 2 = col(0,1,0), generates the vector V2 = col(0,1, 0, 0, 0, 0, 0, 0, 0), while the Kronecker product of V 2 and vi generates = 
col(0, 0, 0, 0,1, 0, 0, 0, 0). These examples coincide precisely with the mapping (n, m) m (n — l)d. 
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tensor product. But above all, what is important to understand what one is doing, and hence I strongly encourage 
the reader to think deeply about this, and play with some examples to interiorize this tricky point. 

Hence, as promised, the matrix representation of the Liouvillian superoperator admits a very simple coding in 
Mat lab: 

L = -li*kron(I, H) + li*kron(H.', I) + 2*r*kron(conj (J), J) - r*kron(I, J'*J) - r*kron(J.'*conj (J), I), ( 14 ) 

where li is the proper way of writing the imaginary unit in Matlab, while conj(z) is how the complex conjugate of 
z looks in Matlab. 

Once we have the Liouvillian superoperator, the next issue concerns finding the steady state of the system. If the 
Hilbert space dimension is not too large, we can try diagonalizing L fully using [V, D] = eig(L). To access the steady 
state we can just find the index of the zero eigenvalue and get the corresponding column of V, or proceed in a more 
elegant and automatic way, by sorting the order in which the eigenvalues appear. In particular, given the vector of 
eigenvalues A = diag(D), we generate a vector y containing the indices of the permutation which we need to apply 
to reorder the eigenvectors in decreasing real part as [x, y] = sort(real(A),‘descend’), where we additionally get the 
vector of sorted real parts x, which we won’t use; once we have y, we can sort the eigensystem as V = V(:,y) and 
A = A(y), and the first column of the sorted V should correspond now to the steady state po = V(:, 1 ), most likely 
requiring the additional normalization po = po/trace(reshape(po, d, d)) to ensure it has unit trace. Now, for larger 
size problems, we will need to use sparse matrices, in which case the simplest way of finding the steady state would 
be as po = eigs(L, I,‘LR’), which might require additional normalization as in the previous line. The density matrix 
of the steady state can then be found as pg = reshape(po, d, d). 

Even though in most cases the previous way of finding po is enough, it is interesting to know how to implement 
the method which we introduced at the end of the first section, which consisted in replacing one equation of Lpo = 0 
by the normalization condition 7(pii + P22 + ••• + pdd) = 7, where 7 > 0 is a parameter which we can choose as 
we wish. This is easily done in Matlab as follows (there are indeed many different ways of doing this, here I just 
pick the one I find simplest and most direct to code). First, we pick the index I of the diagonal element pu whose 
equation we want to replace, with corresponding superindex si = / + (/ — I)d. Then, we just define the matrix Lg = L, 
and replace the corresponding row as Lg(5/,:) = 7*I(:), which simply puts 7 on the entries multiplying the diagonal 
elements of the density matrix. Then we define the vector wq with a single 7 on element si , which indeed corresponds 
to the superspace basis element multiplied by 7, that is, we can simply code it as wq = 7*I(:, si). Once we have 
done this, pg can be found as pg = inv(Lg)*reg, or by asking Matlab to solve the linear system Lgpg = wq as pg = 
linsolve(Lg, rcg), which uses LU factorization. Note that we need to check that det(Lg) 7^ 0 , as otherwise we will 
have degenerate steady states, and the method will fail. 


IV. DEALING WITH COMPOSITE HILBERT SPACES 

Everything we introduced up to now is general, in the sense that it applies to a system with any Hilbert space H. 
Here I want to discuss some special features that appear when the system is a composition of N simpler subsystems 
(two- or three-level systems, harmonic oscillators, etc...), in which case the Hilbert space has the tensor product 
structure H = Hi (8)^2 (8)... 0 Hat- This is the scenario that we usually find in quantum optics, where typical systems 
are composed of atoms (maybe artificial such as superconducting qubits or quantum dots), and/or photonic, phononic, 
or motional modes. 

Given a basis {\nj)}nj=i^2,...,dj of subspace Hj with dimension dj, a basis of the full Hilbert space can be built as 
{|ni) ( 8 ) 1^2) ( 8 ) ... ( 8 ) \nN)}nj=i,2,...,dj = {|^)}n=i,2,...,d5 where d = di x ^2 x ... x is the dimension of the complete 
Hilbert space. Hence, each value of the index n which we were using in the previous sections corresponds now to some 
multi-index nin2...n]sf labeling the basis elements of the Hilbert space of each subsystem. The point is that in many 
cases (for example when wanting to evaluate partial traces or transpositions) it is important to keep track of all the 
indices, and here I want to discuss how to deal with these issues by using efficient computational tools. 

We have already encountered tensor products before when building the superspace, and we saw how to code them 
efficiently using the Kronecker product; we will again use the kron operation to deal with composite Hilbert spaces, 
but with a few subtle points. Let us start with the simplest structure in the composite Hilbert space: let’s represent its 
basis. Consider again the basis {\nj)}nj=i,2,...,dj of subspace Hj, and define the identity matrix of the corresponding 

dimension, I^-^^ = eye(dj), from which we build the vector representation of \nj) as Then, we 

will represent a basis element n of the complete Hilbert space corresponding to some multi-index nin2...nAr, that is, 
1 ^) = 1^1) C) 1^2) ( 8 ) ... 0 as the vector = kron(vi^\ ...,kron(vi^\kron(vi^2\ ))•••)5 where we use again a 

reversed order in the kron operations for future convenience, see the next paragraph. Taking into account that we still 
want to define as a vector with d — 1 zeros and a one at position n, which is the natural self-representation of the 
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basis elements in the complete Hilbert space, the previous definition fixes the relation between n and the multi-index 
nin2...nAr to 


n = m -h (n 2 - l)di -h (ns - l)did 2 + ... + {un - l)did 2 ...dN-i> ( 15 ) 

Consider now a pure state \a) = complete Hilbert space. As explained, it would be useful 

to be able to move between this expression, and the one making explicit reference to the indices of each subspace, 
that is, |a) = ■■■J2tl=lan^n2...n^ |ni) ( 8 ) |n2) 0 ... 0 |nAr). This is very easy to do in Matlab once 

we have made all the previous definitions. In particular, given the (column) vector representation of the state 
a = col(ai,a2, ...,nd) or better a = [ai; a2;...; in Matlab code, we can transform it into a multidimensional array 
as a = reshape(a, di, ^2,..., d^v), from which an^n2...nN simply accessed as an^n2...nN — ^(^i?'^2, n^v); in the 

following we will use the double dot on top of the bold-faced symbol to denote that it is a multidimensional array. 
The multidimensional array can be transformed back to its vector form in the full Hilbert space just vectorizing it as 
a = a(:). Note that the dimensions (di,d2, ...,dAr) used to reshape the vector and the indices (ni,n2, ...,nAr) of the 
multidimensional array, follow the intuitive order that one would assign; this is thanks to using the reversed order in 
the kron operation, and would not be the case if we would have chosen the intuitive one. 

At this point it should be clear that, given substates G V,j}j=i^2,...,N with vec¬ 
tor representation a^-^^ = col(a^'^\..., a^'^.^) or a^-^^ = in Matlab code, the vector repre¬ 
sentation of the tensor product state \a) = ( 8 ) |a( 2 )) ® ... ® |aW) can be obtained in Matlab as a = 

kron(a^^\ ...,kron(a^^\kron(a*^^\ a*^^^))...), with elements = ani Cin2 ^ given by ( 1 ^. 

We see then that dealing with vectors is not so difficult. Now let’s consider operators, which are a bit more tricky. 
Let’s start with a simple generalization of what we did for vectors. Consider an operator O = Ylnm=i Onm|^)(^| in 
the complete Hilbert space. Given its matrix representation O in Matlab (with elements O^m) we would like to be 
able to retrieve the multi-index elements Onin2...nN]mim2...mN defined from 

di d 2 dN 

0= E E - E Onin 2 • • - nN ',Triirn 2 • • - rnN (|ni) ® |n 2 ) ® ... 0 |njv))((«i’i| 0 ("I 2 I 0 ... 0 {mwl)- (16) 

ni,mi = l n 2 ,m 2 = l nN,'mN = ^ 

Similarly to vectors, this can be done by redefining O as a multi-dimensional array O = 
reshape( 0 , di, (^2,..., div, di, (^2,..., d^v), from which we can then get the desired multi-index elements as 
Onin2...nN]mim2...mN = O(ui, 77 - 2 ,..., uat, uii, m2,..., uiAr). Wc cau always go back to the original matrix repre¬ 
sentation in the complete Hilbert space by reshaping the multidimensional array as O = reshape( 0 , d, d). 

Sometimes it is useful to have access to a different set of multi-index elements Onimi;n2m2;...;nA^miv defined by 

di d 2 dN 

0 = E E - E OniTTii ;n 2 m 2 ;...;nArmAr |ni)(mi| 0 |n2)(m2| 0 ... 0 |nAr)(mAr|. ( 17 ) 

ni ,mi=l 712,7712 = 1 777 V, 7 n 7 V = l 

For this, the best is first building the multidimensional array O as we explained above, and then use the extremely 
useful permute operation, which allows to permute indices of multidimensional arrays. In particular, defining another 
multidimensional array O = permute( 0 , [I, A" + 1 , 2 , A + 2 ,..., A, 2 A]), we then access the desired multi-index ele¬ 
ments as Oni7ni;772 7772;...;nAr777Ar = 0 (ni, mi, 77 - 2 , m2,..., utv, m^v). In the following we will use the notation O for the 
multidimensional array corresponding to the order Oni772...77Ar;777 i7772...mAr Ide multi-index elements, and the notation 
O for the one corresponding to the Oni777i;77 2 7772;...;nAr777Ar order. 

Imagine now that we are given a set of operators acting on the subspaces {'Hj}j=i^2,...,N, with 

corresponding matrix representations {Oj}j=i^2,...,N‘ Our goal now is finding the different representations of 
the operator O = 0 i( 8 ) 02 C)...( 8 ) 0 Ar acting on the complete Hilbert space. Our starting point will be the 
Kronecker product O = kron(OAr, ...,kron(03,kron(02, Oi))...); even though this is a d x d matrix containing 
the elements of the representation of O in the complete basis {\n)}n=i^2,...,d^ h is easy to see that these ele¬ 
ments are not ordered in the right way, and hence, in order for O to coincide with the proper matrix repre¬ 
sentation of O we need to reorder its elements. To this aim, we first build the multidimensional array O = 
reshape(0,di,di,d2,d2, ...,dAr,dAr) which has 0 (ni, mi, 77.2, m2,..., ttat, mAr) = 077 i 7 ni; 7727772 ;...; 77 Ar 7 nAr its elements. 
Then, we can reorder its indices as O = permute( 0 , [ 1 , 3 ,..., 2 A — I, 2 , 4 ,..., 2 A]), creating a multidimensional array 
which has 0(77-1,77-2,..., tt-at, mi, m2, ...,mAr) = Otii772...77at; 777i7772...miv il^ elements. Finally, we find the matrix repre¬ 
sentation of O in the complete Hilbert space as O = reshape( 0 , d, d). Hence, essentially we have followed the path 
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of the previous paragraphs, but in reverse order. I hope all the manipulations have served to gain intuition about the 
operations kron, reshape, and permute, which allow for efficient and clean ways of representing Hilbert space objects 
in the computer. 

Finally, I would like to discuss two operations very relevant in the context of composite Hilbert spaces: the partial 
transposition and the partial trace of an operator. Consider an operator O acting on the complete Hilbert space. We 
define the operator corresponding to a partial transpose of O with respect to subspace Hj as 

di dj dN 

= y] ... y] ... y] Oni...mj...nN\mi...nj...mN{\'ni)®-®\nj)®-®\nN)){{mi\®...®{mj\®...®{mN\) 

ni,mi = l nj,mj = l nN,'mN = ^ 
di dj dN 

= E - E - E Om |ni)(mi| (g)... 0 \nj){mj\ 0 ... 0 |nAr)(mjv|. (18) 

ni,mi = l nj,mj = l nN,'fTiN = ^ 

The different multidimensional array representations of this operator are easily obtained in Matlab from the ones of 

••T- •• .... X”. 

the original operator as O ^ = permute( 0 , [1,..., j — 1, + ..., j, +1,..., 2A^]) or O ^ =permute( 0 ,[l,..., 2 (j — 

l),2jf, 2jf — 1,...,2A^]), and its matrix representation in the complete Hilbert space is then obtained as = 

.. T 

reshape (0 \d,d). This is intuitively generalized to the case in which we want to perform partial transposition 
with respect to several subspaces; for example, if we want to transpose subspaces Hj and Hi {j < /), then 


6^^^^ = permute(6, [1,..., j, N j,l,N Hi, j, N H j H 1..., I, N Hi Hi, 2A^]), (19a) 

= permute(0, [1,..., 2(j - 1), 2j, 2j - 1,..., 2(1 - 1), 21, 21 - 1,..., 2N]), (19b) 

qTjTi = 2-eshape(6 ^ \d,d). ( 19 c) 


As for the partial trace over subspace HLj, denoted by tr^jO} or , it is defined as the operator 




ni,mi = l n7v?i^7V = l = ^ 


Eo 


ni mi;... ;nj nj;... ;nA/^ mAT 


ni){mi\ 0 ... 0 \nj-i){mj-i\ 0 \nj+i){mj+i\... 0 |njv)(mjv| 


ni,mi = l n7v,i^7V = l = l 


E^ 


ni...nj ...nN]'m'i---'nj ...rriN 


\nj+i) 0 ... 0 |niv))((TOi| 


that is, as the operators with elements corresponding to the contraction of the indices associated to the Hj subspace. 
We can find again the different representations of this operator efficiently using Matlab’s built-in functions. In 

particular, we can find the multidimensional array O corresponding to this operator as follows (we proceed by 
sequentially updating its definition): first, we bring the indices that we want to contract to the end of the array, 

.. Sj} 

O = permute( 0 , [1,..., j —1, j + 1,..., A^+j —1, A^+j + 1,..., 2 N,j, NHj]); defining the dimension of the total Hilbert 
space after tracing out the desired subspace by d^^^ = di x ... x dj-i x x ... x dw, we reshape the previous array as 

a X d‘j dimensional matrix, 6^"^^ = reshape(6^'^^, d|); given the identity matrix of dimension dj denoted 

by = eye(dj), in terms of the previous matrix the contraction we are looking for is just = 6^' ^^*I(J)(:); the 

previous operation leaves us with a column vector with components, which we can finally reshape to give the 

.. {j} .. {j} 

multidimensional array we are looking for, O = reshape (0 , di,..., dj_i, d^+i,..., dw, di,..., dj_i, d^+i,..., dw); 
finally, we find the matrix representation of the operator in the (remaining) complete Hilbert space as O = 

reshape (0 , d^-^^). Of course, a similar trick can be done starting from the other multidimensional array O; 

also, the method is straightforwardly generalized to when we want to trace out several subspaces at once (maybe it’s 
a good thing to try these two things out as an exercise, to really discover if you understood all these constructions 
properly). 

Let me finally remark once more that the only reason why we have choosen the reversed order in the kron operation 
is to make the coding simpler in Matlab. In particular, by just sticking to the simple rule “every time a tensor product 
appears, it is coded as the kron product in the reversed order”, the rest of manipulations are compactly and intuitively 
coded in Matlab as shown above, what would not be the case otherwise. 
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V. AN EXAMPLE: THREE-LEVEL CASCADE SYSTEM INTERACTING WITH TWO QUANTIZED 
OPTICAL MODES 

In order to fix ideas, let’s consider one example consisting in a three-level cascade S system interacting with two 
driven modes of a cavity which we call a and 6, see Fig. Let us first discuss the structure of this system’s Hilbert 
space as well as a convenient way of writing the master equation governing its evolution, and then we will show how 
to code what we need in Mat lab. 


A. Hilbert space structure and master equation 


The complete Hilbert space of this system can be written as = He. ^ Ha ^ Hh- He is the subspace of the S 
system, with basis {|i)}j=i,2,3, which allows us to define the operators &jk = b)(^|- Ha and Hh are the subspaces of 
the cavity modes, both spanned by Fock states {|n)}^=o,i,2,...: from which we define the basic annihilation operator 
a = l)(^l mode a, and similarly for mode 6, whose corresponding annihilation operator we denote by 

h. Note that Ha and Hh are Hilbert spaces of infinite dimension, but the computer can only deal with finite dimension; 
hence, we need to truncate the Fock state bases to a certain maximum photon number, which we will denote by Na 
and Nh for the corresponding modes, leading to finite-dimensional bases {\n)}a=o,i,2,...,Na {\n)}n=o,i,2,...,Nb which 

approximately span Ha and Hb^ respectively. 

As shown in the figure, we order the states of the S system such that | 3 ) corresponds to the excited state, |2) to 
the middle one, and |1) to the ground one, taking the energy origin in the middle state; we name uji2 and 0023 the 
frequencies of the corresponding transitions. Mode a connects the | 1 ) ^ |2) transition and has resonance frequency 
= ^12 — Aa, detuned by from the transition of the E system. Mode b connects the |2) ^ | 3 ) transition and has 
resonance frequency ujb = 0023 — A5. We assume that both modes are driven by resonant lasers and decay through the 
partially reflecting mirror at rates 7^ and 75 (the other mirror is assumed to have perfect reflectivity, although that’s 
not important for this simple example). Levels | 3 ) and |2) of the E system might decay through modes different than 
the cavity ones, what causes them spontaneous emission to levels |2) and |1), respectively, at rates 723 and 712. All 
these processes are captured by the following master equation in the Schrodinger picture: 


dp 

dt 


-f Ja^a [P] + [p] + 7l2>Ccri2 [P] + l23^a23 [p] ^ 


( 21 ) 


where H{iPj — Hq -H^oupiing -^driving (0 5 with 


Ho = ujaO^a + Wiph + W23CT33 - W12CT11, 

-^coupling = 9 a{a^ ^12 + aa\^) + gb{V'o 20 + bal^), 

HdrWingit) = + e^‘^'’^Sfb), 


(22a) 

(22b) 

(22c) 




FIG. 1 : Sketch of the system used as example: the transitions of a cascade three-level system are coupled to two modes of a 
cavity, as well as to electromagnetic modes outside the cavity which induce spontaneous emission. The cavity modes are driven 
by external resonant lasers and have losses through the partially transmitting mirror. 
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and where we have introduced the notation Cc[p] = 2 c 0 — c^cp — pc^c^ given an operator c. Note that we are not 
writing tensor products explicitly, and hence objects like d(j|2 must be understood as ^{2 we will stick to this 

economic notation except when it can lead to a misunderstanding or we want to show the underlaying tensor product 
structure of the Hilbert space explicitly for some reason. 

Unfortunately, there are two properties of this master equation that make it very difficult to deal with numerically. 
First, it is explicitly time-dependent through ^driving(^)- Second, for large driving amplitudes Ej it is to be expected 
that the cavity modes will get highly populated, and we will not be able to truncate the Fock bases to small enough 
Na and Nij. Both these problems appear typically in many systems, and up to a point can be solved by moving 
to a different picture^, as we will learn now. The idea consists in doing two changes of picture (it can be done at 
once, but it’s more clear in two steps). First, we move to a picture rotating at the laser frequencies (which in this 
example coincide with the cavity frequencies); this is defined by the unitary transformation U{t) = exjp{—iHct) with 
He = + uji)¥b + ujb^ss — SO that the transformed state pu = pU evolves according to^ 

+ laEa [pt/] + IbEb [pu] + 7l2>Ccri2 [M + 723>Ccr23 [M 5 ( 23 ) 

with Hu = Ha ^coupling + where 

HA=Ai,as 3 -Aaau, and = {Saa^ + S:d) + + Sfb). ( 24 ) 


dt 


Hu, Pu 


Hence, we see that in this picture the master equation becomes time-independent. From this new picture we move to 
another one in which, in loose terms, the photons generated by the coherent drivings are already taken into account, 
so that we don’t need to ‘count’ them in the simulation. More specifically, this picture is defined by the unitary 
(displacement) transformation D[a{t), / 3 {t)] = ex.p[a{t)d^ — a*{t)d -h l 3 {t)¥ — ;d*(t)6], which depends on two time- 
dependent amplitudes a{t) and / 3 (t) that will be chosen later. In this case, the transformed state pu = puD 
evolves according to^ 


dpD 

dt 


HD{t),PD 


+ [Pd] + Jb^b [Pd] + 7l2>Ccri2 [Pd] + J23^a23 [Pd] 


{Ea - laOi - a)d) + {Eb - ^bp - - H.C., pu 


( 27 ) 


where Hjjit^ — Ha T -^coupling T with 

^Rabi = 9a[o!*{t)cFi2 + a{t)a\ 2 ] gb[P*{t)d2Z + 


( 28 ) 


® Recall that given the master equation (^, where the Hamiltonian can even be time-dependent, and a general time-dependent unitary 
it is simple to prove that the transformed state pu = pJJ evolves according to the master equation 

= -i Pt/] + - JuJuPu - PuJuJu), (22d) 

with new Hamiltonian Hu = HU + jd£)lJ and jump operators Ju = UHU. 

^ Note that we use here jdt — iHcU^, and 

U^dU = U^bU = e-'^^btb, ai 2 U = and a 23 U = e-'^^^^a 23 , (22e) 

easy to prove from the Baker-Campbell-Haussdorf lemma 

OO -J 

e^ie-s = -IB,[B,...,[B,A]...]], (22f) 

^^ 

n n 

and the commutators [a,a’*’] = 1 and [d-jk^di^] = — ^mjdik- 

® This can be proved by using 

= {ada+a*da* +( 25 ) 

= ^a*a - Aat + /3*S - I3b^ + ^ 

together with 

D^dD = d + a{t) and D^bD = b + p{t). (26) 
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This master equation suggests choosing a and (3 such that its last term is cancelled, that is, as solutions of d = Ea—la^ 
and p = £b- 76/3: 


a{t) = Q;( 0 )e-^“* + — (l - 6“^'*) and ^(t) = / 3 ( 0 )e-^'>* + — (l - . ( 29 ) 

la lb 

Note that with this change of picture we have introduced time dependence in the master equation; however, since 
we are only interested in the long-time behavior of the system (steady state), and moreover, we can choose a and 
(3 at will, we can take the t ^ 7a: 76 limit in the previous equation, in which case the displacements become time 
independent, a = Eajla and (3 = ShMhi and so does the master equation, which takes the final form 

3- 3- ^h^h[pD\ + 721>^cri2 [Pi^] + 732>Ccr23 5 (30) 

with H]J = -^coupling 3 ~ -^Rabi: bciug 


dpD 

dt 


Hd,Pd 


Ha = A^dss - A^dii, ( 31 a) 

-^coupling — 3 ’i 2 3 ~ Qjd’^2) 3 r ^’23 3 ~ ^^23), ( 31 b) 

-^Rabi = (f^a'^12 + ^a^U) + (OJ(723 + ( 31 c) 

where we have introduced the Rabi frequencies fta = Pa^alla and = Pb^bllb- 


B. Coding the problem in Matlab 


In the following we will learn how to code the previous problem in Matlab, with the aim of finding the steady state 
of master equation ( 30 ), and compute certain interesting objects and quantities derived from it. 

One usually starts by defining the basic operators in the complete Hilbert space. For this, we first need to choose 
the bases of the different subspaces and order their elements; in our case, we take the bases that we introduced 
at the beginning of the previous section, ordered as we did (in increasing number of their excitation number). In 
particular, the basis {| 1 ), | 2 ), | 3 )} associated to the energy levels of the S system spans Hs, while the Fock bases 
{| 0 ), | 1 ),..., |A^a)} and {| 0 ), | 1 ),..., span Ha and Hb, respectively. Defining the identity of dimension 3 , the 

representation of eigenvector | 1 ) G Hs corresponds to its first column, while the one of | 3 ) G He to its third column. 
Similarly, defining the identity of dimension + 1 , the representation of Fock state | 0 ) G Ha corresponds to its first 
column, while that of \Na) G Ha to its last column, and the same for mode b. 

As for the basic operators, let’s start from the ones acting on the cascade subspace, the transition operators djk = 
\j){k\. Given the vector representations {vi, V2, V3} of the basis elements {|1), |2), | 3 )}, the matrix representation of 


these operators is obtained in Matlab as ajk 


= 


which has a single one at position {j^k). As for the bosonic 


operators d and 6, note that their matrix elements are {amn = (^|d|n) = i/n6m-\-i,n}n,m=o,i,...,Na {bmn = 
{m\b\n) = where in these expressions the basis vectors are Fock states in the corresponding 

subspaces; hence, their matrix representations are 


/O 1 0 

0 0 V2 

Vo 0 0 


\ 


and b : 


VKJ 


/o 1 0 
0 0 \/2 

Vo 0 0 


0 \ 
0 




( 32 ) 


with the square root of the excitation numbers in the first upper diagonal. Note that these are the representations 
of the operators in their respective subspaces, and they have to be (tensor) multiplied by the identity in the rest 
of subspaces to get their representations in the complete Hilbert space, e.g., G a G in the case of the 

annihilation operator of the a cavity mode. 

With all these considerations and the general constructions of the previous sections, we can start writing the Matlab 
code. In the following, we will go through the main parts of the code, explaining it and writing it explicitly so that it 
can be copied directly to a Matlab script (a simple text file saved with “.m” extension); in any case, the whole script 
can be found as part of the supplemental material. 

We start by giving values to the model parameters A^, A5, ga, Pb^ 712, 723, 7a: 76: ^a: and fib' 
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Deltaa = 0; yodetuning of mode a 

Deltab = 0; yodetuning of mode b 

ga = 1; yocoupling mode a 
gb = 1; yoCoupling mode b 

gammal2 = 1; ^spontaneous decay rate from 2 to 1 
gamma23 = 1; ^spontaneous decay rate from 3 to 2 
gammaa = 3; ^cavity damping rate of mode a 

gammab = 3; ^cavity damping rate of mode b 

Omegaa = 20; “/oRabi frequency driving transition 1-2 
Omegab = 5; “/oRabi frequency driving transition 2-3 

Note that anything written after the symbol (in the same line) is understood as a comment by Matlab, and 
not executed. On the other hand, the semicolons prevent the expression from appearing in the main command 
window (try removing one, and you’ll see how the output of the line is printed on screen). We have chosen simple 
values for the parameters, leading to intuitive physical behavior of the system. In particular, the cavity modes are on 
resonance with their corresponding transitions, the couplings and spontaneous emission rates are on the same order, 
but smaller than the damping through the mirrors, and the Rabi frequencies are the dominant parameters, but with 
the lower transition driven more strongly. Under such conditions, it is to be expected that the population of the S 
system will be almost equally distributed between its ground and middle states, with just a little bit in the excited 
state, and this is exactly what we will see later. 

Let’s now define the parameters related to the dimension of the Hilbert space: 

Na = 4; yoFock basis truncation for mode a 

Nb = 2; yoFock basis truncation for mode b 

dima = Na+1; ^dimension of mode a Hilbert space 

dimb = Nb+1; ^dimension of mode b Hilbert space 

dims = 3; ^dimension of cascade system Hilbert space 

dimtot = 3*dima*dimb; yodimension of the total Hilbert space 

Note that one needs to check that the truncations Na and are enough, by going to larger numbers, and 
confirming that the quantities of interest have converged. 

Now we can start defining the matrix representations of the different operators. We start with the identity operators 
in the different spaces: 

la = speye(dima); yoidentity on mode a subspace 

Ib = speye(dimb); yoidentity on mode b subspace 

Is = speye(3); ^identity on cascade system subspace 

Itot = speye(dimtot); ^identity on the complete Hilbert space 

Note that we have chosen to define them in sparse form to save memory (in full form we would just replace speye 
by eye). The annihilation operators for the cavity modes are then written in their respective subspaces as 

a = spdiagsCsqrt(0:Na)U1jdima,dima); 
b = spdiagsCsqrt(0:Nb)U1jdimbjdimb); 

in sparse form, or 

a = diagCsqrt(1:Na),1); 
b = diagCsqrt(1:Nb),1); 

in full form. In the complete Hilbert space Hs ^ 'Ha ^ Hb^ these operators are coded as 

a = kron(Ib,kron(a,Is)); 
b = kron(b,kron(Ia,Is)) ; 

as we learned in Section lYl 

In order to code the transition operators ajk of the cascade system, it is convenient to first define the vector 
representation of its basis elements, what we do as 
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vl = Is(:,l) 
v2 = Is(:,2) 
v3 = Is(:,3) 


“/oground state of the cascade system 
“/omiddle state of the cascade system 
“/oexcited state of the cascade system 


Once we have the basis vectors, we can code the transition operators in the complete Hilbert space as 

sll = kron(Ib,kron(Ia, vl*vl O ) ; yoSigma_{ll} 
s22 = kron(Ib,kron(Ia, v2*v2O ) ; yoSigma_{22} 
s33 = kron(Ib,kron(Ia,v3*v30); yoSigma_{33} 
sl2 = kron(Ib,kron(Ia,vl*v20); yoSigma_{l2} 
sl3 = kron(Ib,kron(Ia, vl*v30 ) ; yoSigma_{l3} 
s23 = kron(Ib,kron(Ia, v2*v30 ) ; yoSigma_{23} 

Having the matrix representations of the fundamental operators, we are in conditions to code the Liouvillian as a 
matrix in superspace. For this, it is convenient to first code the Hamiltonian, what we do as 


yoBuild the term containing the detunings: 

HDelta = Deltab*s33-Deltaa*sll; 
ythe coupling terms: 

Hcoupling = ga*(a^*sl2+a*sl2O + gb*(b^*s23+b*s230; 
yand the Rabi terms: 

HRabi = (conj(Omegaa)*sl2+0megaa*sl2O + (conj(Omegab)*s23+0megab*s23O; 
yfrom which we build up the total Hamiltonian: 

H = HDelta+Hcoupling+HRabi; 


Next we code the different dissipative pieces of the Liouvillian. As we learned in Section HL this can be done as 


yoDamping term of mode a: 

La = gammaa*(2*kron(conj(a),a)-kron(Itot,a^*a)-kron(a.^*conj(a),Itot)); 
ydamping term of mode b: 

Lb = gammab*(2*kron(conj(b),b)-kron(Itot,b^ *b)-kron(b.^ *conj(b),Itot)); 
yradiative decay of the lower transition of the cascade system: 

L12 = gammal2*(2*kron(conj(sl2),sl2)-kron(Itot,sl2 ^ *sl2)-kron(sl2.^ *conj(sl2),Itot)); 
yradiative decay of the upper transition of the cascade system: 

L23 = gamma23* (2*kron(conj (s23) , s23)-kron(Itot, s23^ *s23)-kron(s23.’*conj (s23) ,Itot)); 


Once we have the Hamiltonian and the dissipative pieces, we then build the total Liouvillian as 


L = -li*kron(Itot,H) + li*kron(H. \ Itot)+La+Lb+L12+L23; ytotal Liouvillian 


as given by expression (14). 


At this point we have managed to code the matrix representation of the Liouvillian in superspace. Now, we proceed 
to evaluate its steady state in the different ways that we introduced in Section [Hi] . For each method, given the steady 
state which we denote here by ps, we compute the populations tr{djjps}, trjd^dps}, and tr{6'*'6ps}. At the end we 
will see that all the methods give the same populations. 

As a first method we find the eigenvector with zero eigenvalue via sparse diagonalization, as explained in Section 
imi The code looks like 


[rhoSl,lambdaO] = eigs(L,1,^LRO; yfind eigenvector with largest real part 
eigenO = lambdaO ycheck that the eigenvalue is 0 

rhoSl = reshape(rhoSljdimtotjdimtot); yreshape eigenvector into a matrix 
rhoSl = rhoSl/trace(rhoSl); ynormalize 

Popl = [trace(sll*rhoSl) trace(s22*rhoSl) trace(s33*rhoSl)... 

traceCa^*a*rhoSl) trace(b^*b*rhoSl)]; yevaluate populations 

The second line prints out the eigenvalue of the Liouvillian matrix with the largest real part, which should appear in 
Matlab’s command window as 


eigenO = 

-9.6655e-15 - 7.7851e-15i 
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Note that this is basically zero within the numerical error, just as expected. Note also that we have introduced the 
three dots which is just a way of telling Matlab that the expression is too long, and it continues in the next 

line, so lines connected by three dots are understood as a single line by Matlab. 

Let’s consider now the method which uses the full diagonalization of the Liouvillian matrix. We can code it as 

tic “/o St art counting time 

[V,D] = eigCfull(L)) ; “/ofind full eigensystem of the Liouvillian 
t_FullDiag = toe “/otime lapsed since the previous tic 
lambdav = diag(D) ; “/oeigenvalues 

“/oSort eigenvalues in descending order of the real part: 

[x,y] = sort (real(lambdav) , ^descendO ; “/oy stores the permutation to rearrange 

V = V(:,y); “/oSort the eigenvectors 

lambdav = lambdav(y); “/osort the eigenvalues 

eigenv = lambdav(l:5) “/oShow the first 5 eigenvalues 

rhoS2 = V(:,l); “/othe steady state should be the first eigenvector 

rhoS2 = reshape (rhoS2 jdimtot jdimtot) ; “/oreshape it as a matrix 

rhoS2 = rhoS2/trace(rhoS2) ; “/onormalize it 

Pop2 = [trace(sll*rhoS2) trace(s22*rhoS2) trace(s33*rhoS2)... 

trace(a^*a*rhoS2) trace(b^*b*rhoS2)] ; “/oevaluate populations 

We have introduced the functions tic and toe, which allow to check the time that Matlab needed to evaluate the 
instructions between them. The rest just follows the recipe that we learned in Section in This piece of the code 
prints out the following lines in Matlab’s command window: 

t_FullDiag = 

35.713 

eigenv = 

1.1758e-15 - 1.9065e-14i 
-1.0631 + 3.1308e-14i 
-1.5594 - 20.62i 
-1.5594 + 20.62i 
-1.5596 - 20.617i 


The first quantity is the time needed to perform the full diagonalization of the Liouvillian (in seconds); you can 
check when running the whole code that 35 seconds is approximately 90% of the whole time. The next quatities 
correspond to the eigenvalues with the largest real part; note that only one is zero (within the numerical error), and 
the rest have all negative real parts, so we see that we really have a unique steady state. You can check that the 
instruction eigs(L,5, MRO gives the same 5 eigenvalues, but 60 times faster, showing the power of working with 
sparse matrices. 


As a final method, we code the one in which one equation defining the steady state is substituted by the normalization 

It can be done as follows: 


condition, as explained in Section III 


Isuper = eye(dimtot*dimtot) ; “/odefine the identity in superspace 

1 = 1; “/opick the index of the diagonal element whose equation we want to replace 

si = 1+(1-1) *dimtot; yoCorresponding index in superspace 

gamma = 1; yoConstant by which we multiply the normalization condition 

LO = full(L); “/oWe first define LO as the Liouvillian in non-sparse form 
“/oAnd then replace the chosen row by the part of the normalization condition: 
L0(sl,:) = gamma*Itot(:); 

yoDefine the vector encoding the other part of the normalization condition: 
wO = gamma*Isuper(:,si) ; 

rhoS3 = L0\w0; “/oSteady state in terms of the inverse of LO 

rhoS3 = reshape(rhoS3,dimtot,dimtot); ^reshape it as a matrix 

tr3 = trace(rhoS3) ycheck the trace, which should be 1 by construction 

rhoS4 = linsolve(L0,w0); ^steady state using the linear solver of Matlab 
rhoS4 = reshape(rhoS4,dimtot,dimtot); ^reshape it as a matrix 
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tr4 = trace(rhoS4) yoCheck the trace, which should be 1 by construction 

yWe finally evaluate the populations with both states 

Pop3 = [trace(sll*rhoS3) trace(s22*rhoS3) trace(s33*rhoS3)... 

traceCa^*a*rhoS3) trace(b^*b*rhoS3)]; ^populations from rhoS3 
Pop4 = [trace(sll*rhoS4) trace(s22*rhoS4) trace(s33*rhoS4)... 

trace(a^*a*rhoS4) trace(b^*b*rhoS4)]; ^populations from rhoS4 


Note that we find the steady state by solving its defining equation in the two different ways explained in Section [III 
either by inversion of the modified Liouvillian or using Matlab’s linear solver. The code prints out in Matlab’s 
command window the trace of the density matrices obtained through both methods, which should be 1 by 
construction. You can check that this is indeed the case. 


Next in the code, we evaluate some reduced states as an example of how to code the partial trace. We start from 
the steady state evaluated via sparse diagonalization, rearranged as a multidimensional array as 

rhoMDA = reshape(rhoSl,dims,dima,dimb,dims,dima,dimb); 

From this, we find the reduced state of the cavity modes by tracing out the E system as 

rhoab = permute(rhoMDA, [2,3,5,6,1,4] ) ; yomove cascade indices to the end 

rhoab = reshape(rhoab,dima*dimb*dima*dimb,dims*dims); ^reshape as a matrix 

rhoab = rhoab*Is(:); ytrace out the cascade subspace 

“/oReshape the superspace vector as a multidimensional array: 

rhoab = reshape(rhoab,dima,dimb,dima,dimb); 

yand build the reduced density matrix in the a+b subspace: 

rhoab = reshape(rhoab,dima*dimb,dima*dimb); 

We can also trace out the cavity modes, to find the reduced state of the E system: 

rhos = permute(rhoMDA,[1,4, 2,3,5,6]); ymove cavity indices to the end 
rhos = reshape(rhos,dims*dims,dima*dimb*dima*dimb); ^reshape as a matrix 
lab = eye(dima*dimb); ^define identity matrix in the a+b subspace 
rhos = rhos*Iab(:); “/otrace out the cavity modes 

“/oReshape superspace vector into the reduced matrix in the cascade subspace: 
rhos = reshape(rhos,dims,dims); 

Finally, we find the reduced state of each cavity mode from their combined reduced state found before, first for 
mode a: 


“/oReshape their combined state as a multidimensional array: 
rhoa = reshape(rhoab,dima,dimb,dima,dimb); 

rhoa = permute(rhoa,[1,3,2,4]); ymove indices of the b subspace to the end 
rhoa = reshape(rhoa,dima*dima,dimb*dimb); ^reshape as a matrix 
rhoa = rhoa*Ib(:); ytrace out the b mode 

“/oReshape the superspace vector into the reduced matrix in the a subspace 
rhoa = reshape(rhoa,dima,dima); 

and then for mode b: 

rhob = reshape(rhoab,dima,dimb,dima,dimb); 

rhob = permute(rhob, [2,4,1,3]); ymove indices of the a subspace to the end 
rhob = reshape(rhob,dimb*dimb,dima*dima); ^reshape as a matrix 
rhob = rhob*Ia(:); ytrace out the a mode 

“/oReshape the superspace vector into the reduced matrix in the b subspace: 
rhob = reshape(rhob,dimb,dimb); 

Now that we have found the reduced steady states, let’s compute the populations from them. 
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yoDefine operators in their respective subspaces: 

ar = diagCsqrt(1:Na),1); ^annihilation operator in the a subspace 
br = diagCsqrt(1:Nb),1); ^annihilation operator in the b subspace 
sllr = vl*vl^; yoSigma_{ll} 
s22r = v2*v2^; yoSigma_{22} 
s33r = v3*v3^; yoSigma_{33} 

PopReduced = [trace(sllr*rhos) trace(s22r*rhos) trace(s33r*rhos)... 

traceCar^*ar*rhoa) trace(br^*br*rhob)]; ^populations 

Then, we build a matrix containing the populations from all the methods as columns (note that we take the real 
parts, so that the imaginary parts are not printed out to save space on screen, but check yourself that the latter are 
zero as they should be): 


Pop = real([Popl; Pop2; Pop3; Pop4; PopReduced]O 
which printed in Matlab’s command window reads: 


Pop 


0.45882 

0.48438 

0.056796 

0.019165 

0.0012705 


0.45882 

0.48438 

0.056796 

0.019165 

0.0012705 


0.45882 

0.48438 

0.056796 

0.019165 

0.0012705 


0.45882 

0.48438 

0.056796 

0.019165 

0.0012705 


0.45882 

0.48438 

0.056796 

0.019165 

0.0012705 


showing that all the steady states give exactly the same populations. Note that the populations of the S system are 
what we were expecting from the system parameters. On the other hand, note also that we have computed is not 
the true cavity populations, since our state is not in the Schrodinger picture, but in a displaced picture where the 
external driving is subtracted. Taking into account that the steady state is in the Schrodinger picture, 

we can get the true cavity populations as 


aUDp^D'^lj'^} = d) alJDp^} = tr{(d^ + a*){a + (a)ps} = \a\^ + trjd^dps} + 2 Re{(a*tr{dps}}, (33a) 

Dp^} = tr{(6^ + a'){h + (a)ps} = \ol\^ + tT{Pbps} + 2 Re{(a*tr{6ps}}, (33b) 


which we compute in Matlab as 

alpha = Omegaa/ga; 
beta = Omegab/gb; 

Popa = Popl(4)+conj(alpha)*alpha+2*real(conj(alpha)*trace(a*rhoSl)) 
Popb = Popl(5)+conj(beta)*beta+2*real(conj(beta)*trace(b*rhoSl)) 

printing out the following result in Matlab’s command window: 

Popa = 

399.66 + 3.9078e-16i 
Popb = 

24.961 + 2.1115e-16i 


Hence, we see that with such strong drivings, the S system doesn’t change too much the cavity populations from 
their values expected in the absence of coupling, = 400 for mode a and = 25 for mode b. 

Finally in the code, we proceed to check the entanglement between various bipartitions of the complete system, 
what will give us a perfect excuse to compute some partial transpositions. Given the state p of a system whose Hilbert 
space we divide in two as Ha G Hb, a necessary condition for it to be separable with respect to that bipartition is 
that the partial transpose p^^ is semi-positive definite, that is, it has only positive or zero eigenvalues. Given the 
eigenvalues {A^jn of we can evaluate the level of violation of such condition via the logarithmic negativity 
F^ln = log[l + “ ^n)]? which is one of the most common entanglement measures available for mixed states. 

In the following we evaluate this quantity for various bipartitions of our system. 

Let’s start with the entanglement between the S system and the cavity modes. We can find the corresponding 
logarithmic negativity as 
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yoGiven the full state as a multidimensional array, 
ywe first transpose the cascade subspace: 
rhoT = permute(rhoMDA,[4,2,3,1,5,6]); 

rhoT = reshape(rhoT,dimtot,dimtot); ^reshape it as a matrix 
Teigenv = eig(rhoT); ^compute its eigenvalues 

logNeg = log(l+sum(abs(Teigenv)-Teigenv)); ^compute the log negativity 

Let’s compute now the entanglement between the cavity modes, what we do as 

yFirst reshape the reduced state of the cavity modes 

yas a multidimensional array: 

rhoabT = reshape(rhoab,dima,dimb,dima,dimb); 

rhoabT = permute(rhoabT,[3,2,1,4]); ^transpose the a mode subspace 
rhoabT = reshape(rhoabT,dima*dimb,dima*dimb); ^reshape it as a matrix 
abTeigenv = eig(rhoabT); ^evaluate its eigenvalues 

logNegab = log(l+sum(abs(abTeigenv)-abTeigenv)); ^compute the log negativity 

Finally we check the entanglement between the S system and each of the cavity modes individually. We start with 
the a mode as 

yFirst we need the reduced state of the cascade system and the a mode. 

^Starting from the complete state as a multidimensional array, 
ywe move the indices of the b mode to the end: 
rhoas = permute(rhoMDA,[1,2,4,5,3,6]); 

rhoas = reshape(rhoas,dima*dims*dima*dims,dimb*dimb); ^reshape as a matrix 

rhoas = rhoas*Ib(:); ytrace out the b subspace 

yand reshape the superspace vector as a multidimensional array: 

rhoas = reshape(rhoas,dims,dima,dims,dima); 

rhoasT = permute(rhoas,[3,2,1,4]); ^transpose cascade subspace 
rhoasT = reshape(rhoasT,dima*dims,dima*dims); ^reshape it as a matrix 
asTeigenv = eig(rhoasT); ind eigenvalues 

logNegas = log(l+sum(abs(asTeigenv)-asTeigenv)); ^compute the log negativity 
and then for the b mode as 

yFirst we need the reduced state of the cascade system and the b mode. 

^Starting from the complete state as a multidimensional array, 
ywe move the indices of the a mode to the end: 
rhobs = permute(rhoMDA,[1,3,4,6,2,5]); 

rhobs = reshape(rhobs,dimb*dims*dimb*dims,dima*dima); ^reshape as a matrix 

rhobs = rhobs*Ia(:); ytrace out the a subspace 

yand reshape the superspace vector as a multidimensional array: 

rhobs = reshape(rhobs,dims,dimb,dims,dimb); 

rhobsT = permute(rhobs,[3,2,1,4]); ^transpose cascade subspace 
rhobsT = reshape(rhobsT,dimb*dims,dimb*dims); ^reshape it as a matrix 
bsTeigenv = eig(rhobsT); ind eigenvalues 

logNegbs = log(l+sum(abs(bsTeigenv)-bsTeigenv)); yocompute the log negativity 

Collecting all the logarithmic negativities in a single vector as 

LogNegativities = [logNeg; logNegab; logNegas; logNegbs] 

we get the following printed out in Matlab’s command window: 

LogNegativities = 

0.0025892 - 1.061e-16i 
2.027e-07 - 9.707e-17i 
0.0017957 - 1.2758e-16i 
9.2002e-05 - 1.5018e-17i 
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This shows that there is indeed entanglement between all the bipartitions, although it is not very big in any case 
(take a bell state as an example, which has logarithmic negativity equal to log 2 ^ 0.7), consistent with the fact that 
the cavity populations are not very much affected by the coupling to the S system. Note in particular that the 
largest entanglement is between the S system and the cavity modes, while there is almost no entanglement between 
the cavity modes themselves. On the other hand, the entanglement of the S system with the a mode is much larger 
than that with the b mode. 
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